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1.  INTRODUCTION 

The  problem  of  selecting  a  suitable  probability  model  for  the 
number  of  cycles  to  crack  initiation  was  discussed  by  Soyer  (1982) , 
(1983) .  There  it  was  reported  that  the  Weibull  and  the  lognormal  dis¬ 
tributions  were  appropriate  models  for  describing  the  distribution  of 
the  number  of  cycles  to  crack  initiation. 

In  this  report,  a  related  problem,  namely,  the  estimation  of  the 
threshold  parameter  in  these  life  models,  is  considered.  The  threshold 
parameter  is  the  minimum  value  that  a  random  variable  can  take.  For 
fatigue  crack  problems  it  is  the  largest  number  of  cycles  before  which 
crack  initiation  is  impossible. 

In  this  paper,  we  discuss  a  Bayesian  approach  for  estimating  the 


threshold  parameter  of  the  Weibull  and  the  lognormal  models. 
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2.  ESTIMATION  OF  THE  THRESHOLD  PARAMETER 
OF  A  WEIBULL  DISTRIBUTION 


The  three-parameter  Weibull  distribution  has  the  density  function 


f(x) 


(J),6,0  >  0,x  >  4)  , 


(2.1) 


where  6  and  0  are  the  shape  and  the  scale  parameters,  respectively,  and  4> 
is  the  threshold  parameter. 

In  this  report,  we  shall,  following  the  Bayesian  paradigm,  assume 
that  prior  information  concerning  the  number  of  cycles  to  crack  initia¬ 
tion  exists,  and  that  such  information  can  be  translated  into  prior  dis¬ 
tributions  for  the  parameters  of  (2.1). 

Let  X  be  a  random  variable  having  the  density  (2.1).  The  dis¬ 
tribution  of  X  is  dependent  on  the  unknown  parameters  4),  S,  and  0.  If 
h(4),3,6)  denotes  the  joint  prior  distribution  of  the  parameters,  then 
given  n  lifetimes  x  =  (x^,...,x^)  and  Jl(4>,S,0;  x)  the  likelihood  func¬ 
tion,  by  Bayes'  theorem  the  posterior  distribution  g(4>,B,0  |  x)  is  pro¬ 
portional  to  h(4),6,0)  •  il(4),6>6;  x) ,  the  constant  of  proportionality 
being  independent  of  the  parameters.  If  3  is  an  element  of  the  param¬ 
eter  space  8  and  0  an  element  of  the  space  (h) ,  then  the  posterior 
distribution  of  <p  can  be  obtained  as 


g]^(4'l5) 


g(4).6,0|x)  dQdS  . 


(2.2) 


The  modal  value  of  the  above  posterior  distribution  can  be  used  as  a 
Bayes  estimate  of  the  threshold  parameter  <p.  This  value  is  also  known 
as  the  generalized  maximum  likelihood  estimator,  and  is  viewed  as  a 
Bayes  estimator  when  a  particular  loss  function  is  not  specified. 


The  aforecited  analysis  poses  several  difficulties,  because  there 
exists  no  closed  form  solution  for  the  integration  in  (2.2);  thus  numeri¬ 
cal  Integration  techniques  will  have  to  be  used. 

We  will  consider  two  different  cases.  In  the  first  case,  inde¬ 
pendent  prior  distributions  will  be  assumed  for  all  three  parameters. 

In  the  second  case,  we  will  consider  a  joint  prior  distribution  for  cp 
and  9,  and  a  prior  distribution  for  S,  which  is  independent  of  the  joint 
distribution  of  (p  and  6. 


2.1  Independent  Priors  on  All  Parameters 

Assume  that  a  sample  of  n  lifetimes  is  observed  from  a  Weibull 
distribution  having  the  density  (2.1).  The  likelihood  function  is 

,6 


£(4i,0,0;  x)  = 


n  j 
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i=l 
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6  J 


(2.3) 


We  assume  a  uniform  prior  over  (0,M)  for  the  threshold  parameter  (});  thus 
the  probability  density  function  for  (j)  is: 

hj^(<P)  =  I  ,  0  <  (f.  <  M  . 

Similarly,  the  prior  on  6  is  assumed  to  be  uniform  over  (Mj^,M2),  giving 


^ 


Mj  <  6  «  , 


as  its  probability  density  function. 

For  the  scale  parameter  0,  we  consider  a  four  parameter  beta  den¬ 


sity,  so  that 
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3'  -  r(g^)r(g2)  g,-Hg2-i 

(b-a) 


(2.3.1) 


where  a  <  9  <  b,  g^,g2  >  0. 

Because  of  independence,  the  joint  prior  distribution  of  ip,  g, 
and  9  can  be  written  as  the  product  of  the  marginal  densities, 

h(<t),g,9)  =  h^((J))h2(6)h3(9)  .  (2 


9  is : 


By  using  Bayes'  theorem  the  joint  posterior  density  of  p,  g,  and 


g((t),g,9|x) 


/^in  n  X.- 
cc  ^  n 

9  .  ,  9 

1=1 


n  x.-d) 

e-ll  -  I  -V- 


gi~l  go~^ 

(9-a)  (b-0)  .  (2.5) 


In  order  to  obtain  the  posterior  distribution  of  the  double 


integral  with  respect  to  g  and  9  given  below  should  be  evaluated: 


82^('t>l21^ 


f/r  iir 


n  X 


(2.6) 


g  -1  g  -1 

(e-a)  (b-0)  d9dg 


a 


There  exists  no  closed  form  solution  for  the  above  integral;  thus 
a  numerical  technique  is  used  to  evaluate  it. 

This  approach  is  applied  to  some  engine  disk  data  presented  in 
the  appendix.  A  uniform  prior  density  is  assumed  for  the  threshold  pa¬ 


rameter  p,  over  (0,700);  this  density  is  indicated  in  Figure  1.  The 
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prior  for  B  is  assumed  to  be  a  uniform  over  the  interval  (1,6),  and  a 
beta  prior  density  over  the  interval  (200,900)  is  considered  for  6, 
with  =  3  and  =  2;  see  Equation  (2.3.1).  The  integral  in  (2.6)  is 

solved  by  using  Simpson's  rule  and  the  posterior  distribution  of  <j)  ob¬ 
tained  from  (2.6)  is  shown  in  Figure  2.  Note  that  although  we  started 
with  a  uniform  prior  on  (j),  the  data  have  enabled  us  to  revise  our  prior 
evaluation,  to  yield  a  posterior  distribution  which  has  a  mode  at  375 
equivalent  flight  hours  (EFH) .  This  value  may  be  used  as  the  Bayes 
estimate  of  the  threshold  parameter. 

The  assumption  of  independent  priors  on  all  the  parameters  may 
not  be  realistic;  we  drop  this  assumption  in  the  next  section. 

2.2  A  Joint  Prior  on  the  ThresholJ 
and  the  Scale  Parameters 

In  this  section  a  joint  prior  distribution  is  considered  for  the 
parameters  0  and  0,  and  an  independent  prior  density  is  considered  for  3. 

The  joint  prior  distribution  considered  here  was  first  suggested 
by  Varde  (1969)  for  the  two  parameter  exponential  distribution.  This 
joint  density  has  the  form: 

h2(<^.e)  cc  •  6(n-<l))  ,  (2.7) 

^  eVTj. 

where 

Cl,  u  >  0 

6(u)  =  ' 

(O,  u  <  0  ; 

0  >  0  ,  and  6  >  0  . 

The  parameters  v  and  x  are  positive  integers  such  that 
V  <  X  <  p/q  ;  y  >  0  ,  y  and  q  are  such  that  q  >  0  .  Note  that  when 


X  =  y  =  0  and  rj  ^  (2.7)  is  equivalent  to  assuming  independent  priors 

on  <p  and  0,  with  the  prior  of  (p  being  diffuse  and  the  prior  of  9  being 
proportional  to  l/o'^'*'^. 

A  plot  of  this  joint  prior  density  for  v  =  1,  q  =  700,  y  =  6500, 
and  X  =  3  is  shown  in  Figure  3. 

If  we  integrate  (2.7)  with  respect  to  0,  the  marginal  prior  for 
the  threshold  parameter  is  obtained  as 

h  (4))  <r  ^  6(n-4))  .  (2.8) 

(y-X(())'^ 

A  plot  of  this  marginal  prior  density  is  shown  in  Figure  4.  As 
can  be  seen  from  this  figure,  the  prior  density  increases  with  (p.  The 
parameter  X  controls  the  slope  of  the  density  function.  For  small  values 
of  X  the  prior  becomes  flatter. 

The  parameter  q  determines  the  range  of  values  of  p  over  which 
its  density  is  nonzero;  that  is,  n  is  the  largest  value  that  p  can  take. 

The  parameters  y  and  v  are  harder  to  interpret,  but  once  the 
parameters  r|  and  X  are  determined,  the  inequality  v  <  X  <  y/n  enables  us 
to  specify  the  values  of  these  parameters. 

An  independent  uniform  prior  density  on  (Mj^,M2)  is  assumed  for 
the  shape  parameter  B,  and  thus  the  joint  prior  on  (}),  B,  and  0  can  be 
written  as: 

h(<{),B,0)  -  ^-(u-W/d  _  (2.9) 

0V-ri  ^2 

By  Bayes'  theorem  the  joint  posterior  density  of  the  three  param¬ 
eters  is  obtained  as: 


Density 
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(2.10) 


To  obtain  the  marginal  posterior  distribution  of  ,  we  need  to 
integrate  (2.10)  with  respect  to  6  and  0.  Here  again,  there  is  no 
closed  form  solution  for  the  posterior  distribution  of  v,  and  therefore 
the  above  integral  is  to  be  solved  numerically. 

This  approach  is  applied  to  the  engine  disk  data  by  using  differ¬ 
ent  forms  of  the  joint  prior  (2.7). 

The  prior  density  for  6  is  assumed  to  be  uniform  over  the  interval 
(1,6).  The  joint  prior  on  and  G  will  have  the  parameters  n  =  700, 

U  =  6500,  V  =  1,  and  A  =  5  (see  Figure  3).  The  marginal  prior  density 
for  is  shown  in  Figure  5.  The  posterior  distribution  of  $  is  obtained 
by  solving  the  double  integral  in  (2.10)  numerically  with  respect  to  6  and 
and  0.  The  posterior  distribution  of  i+i  is  presented  in  Figure  6.  The 
modal  value  is  at  500  EFH. 

Another  joint  prior  considered  for  0  and  0  has  parameters  n  =  700, 
y  =  6500,  v  =  3,  and  A  =  5.  The  prior  on  g  is  again  a  uniform  over  (1,6). 
By  using  numerical  integration  in  (2.10)  the  posterior  distribution  of  <fi 
is  obtained.  The  mode  of  the  posterior  distribution  occurs  at  600  EFH. 
This  posterior  distribution  is  shown  in  Figure  7.  Another  joint  prior 
for  (|)  and  0  has  parameters  n  ~  700,  y  =  6500,  v  =  2,  and  A  =  5.  In  this 
case  the  mode  of  the  posterior  occurs  at  550  EFH  (see  Figure  8) .  Thus 


DCMSITY 


Figure  5.  Marg.inal  prior  density  for  threshold  parameter 
with  n  =  700,  u  =  6500,  v  =  1,  and  X  =  5. 


Marginal  posterior  density  for  threshold 
parameter  (p  assuming  the  joinr  prior  of 
Figure  3. 
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Figure  7.  Marginal  posterior  density  for  threshold 
parameter  <{)  assuming  a  joint  prior  for  (j) 
and  0  with  n  =  700,  p  =  6500,  v  =  3, 
and  A  =  5. 


DENSITY 


it  appears  that  increasing  the  value  of  v  shifts  the  mode  of  the 
posterior  distribution  closer  to  n* 


In  order  to  see  the  effect  of  a  decrease  in  the  value  of  the 
parameter  A,  we  consider  the  joint  prior  on  4>  and  9  with  parameters 
n  =  700,  p  =  6500,  V  =  2,  and  A  =  3.  The  posterior  distribution  of  (p 
is  shown  in  Figure  9.  The  posterior  mode  occurs  at  300  EFH.  That  is,  a 
decrease  in  A,  with  the  other  parameters  held  fixed,  shifts  the  posterior 
mode  closer  to  2ero.  This  can  be  observed  by  comparing  the  modes  in 
Figures  8  and  9. 

Thus  it  appears  that  the  posterior  mode  is  sensitive  to  the 
values  of  the  prior  parameters.  As  we  have  observed,  the  slope  of  the 
marginal  prior  density  of  ({>  can  be  controlled  by  parameter  A.  Small 
values  of  A  give  a  flat  prior  density  for  <p,  whereas  large  values  of  A 
give  more  weight  to  the  higher  threshold  values. 


3.  ESTIMATION  OF  THE  THRESHOLD  PARAMETER 
OF  A  LOGNORMAL  DISTRIBUTION 


The  three-parameter  lognormal  distribution  has  the  density  function 


f(x)  = 


i  ^n(x-(j))- 

2l  a 


(3.1) 


/2'rT  a(x-(p) 

where  cj)  is  the  threshold  parameter,  y  is  the  mean  of  the  random  variable 
2 

Y  =  Zn(x-(p) ,  and  a  is  the  variance  of  Y  .  Given  n  lifetimes  x  = 

(x, ,...,x  ),  the  likelihood  function  is: 

1  n  ’ 


Figure  9.  Marginal  posterior  density  for  threshold 
parameter  <|)  assuming  a  joint  prior  for  4) 
and  6  with  r|  =  700,  u  =  6500,  v  =  2, 
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«.(({), p, a  ;x)  = 


(2Tra2)"/2 


"  1 

n  7  ■  vr 

i=i 


n  |£n(x^-({))-p 
“  T  1 


i=l 


.  (3.2) 


Some  features  of  this  likelihood  function  are  discussed  in  Hill 

(1963).  It  is  shown  that  there  exist  paths  along  which  the  likelihood 

2 

function,  for  any  sample,  tends  to  +»  as  the  triple  (4),p,a  )  approaches 
(X(l),  >  where  is  the  first-order  statistic.  Thus  maximum 

likelihood  estimation,  if  used  here,  would  lead  to  difficulties.  From 
a  Bayesian  point  of  view,  there  is  no  real  problem,  since  the  maximum 
of  the  likelihood  function  is  not  of  interest. 

The  goal  is  to  obtain  the  posterior  distribution  of  the  threshold 
parameter  (p,  and  use  the  mode  of  this  distribution  as  the  Bayes  estimate 
of  d). 


Two  different  cases  are  considered  in  our  analysis.  In  the  first 

case,  independent  prior  densities  are  assumed  for  all  three  parameters. 

2 

In  the  second  case,  a  joint  prior  density  is  assumed  for  p  and  0  , 

whereas  a  prior  density  which  is  independent  of  the  joint  distribution 
2 

of  p  and  a  is  assumed  for  (j). 


3.1  Independent  Priors  on  All  Parameters 

2 

A  uniform  prior  density  over  (0,M)  is  assumed  for  a  .  For  the 

time  being,  assume  independent  priors  for  <|)  and  p.  Let  hj^(d))  and  h.,(p) 

denote  the  prior  densities  for  p  and  p,  respectively.  The  likelihood 

function  is  given  in  (3.2).  By  using  Bayes'  theorem,  the  joint  posterior 

2 

density  of  p,  p,  and  a  is  obtained  as: 
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g((t),u,a  lx) 


«  hj|^(({))h2(u) 


I'J  J 


n/2  I  n 


'A  (v^> 


I  [inCx 

2o  1=1  -v 

e  .  (3.3) 


To  obtain  the  posterior  distribution  of  ij),  (3.3)  is  integrated 

2  2 
with  respect  to  O  and  p.  The  integral  of  (3.3)  with  respect  to  a  can 

be  obtained  in  closed  form.  The  integral 


^  .^>n/2  -ZJ  .n^n(x.-(!.)-M]- 

J  b 


2a  i=l 


do 


(3.4) 


can  be  evaluated  analytically  for  large  M  .  For  large  values  of  M 
(3.4)  becomes: 


--  1 
2  ^ 


n 

■X  I  [iln(x  -4))-p]‘ 
1=1  ^ 


J  -  1|  • 


(3.4.1) 


To  obtain  the  posterior  distribution  of  0,  the  following  integral 
needs  to  be  evaluated: 


/"i' 


®  I  h, (0)h2(u) 


i=l 


--  1 
2  ^ 


[  lX.n(x  -0)-ij]‘ 
vi=l 


du  .  (3.5) 


The?fe  exists  no  closed  form  solution  for  (3.5),  therefore  it  is 
evaluated  by  numerical  integration. 

This  approach  is  applied  to  the  engine  disk  data.  A  uniform 
prior  density  is  assumed  for  p  over  (-100,100).  The  prior  for  0  is  a 
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beta  density  over  (50,750)  with  parameters  =  2  and  %2  ~  ^  Equa¬ 

tion  (2.3.1)].  The  prior  density  of  <])  is  shown  in  Figure  10.  The  pos¬ 
terior  distribution  (3.5)  is  obtained  by  numerical  integration.  The 
mode  of  the  posterior  distribution  is  at  550  EFH.  The  posterior  distri¬ 
bution  of  (J)  is  shown  in  Figure  11. 


2 

3.2  Joint  Prior  on  y  and  a 

2 

The  joint  prior  density  which  is  considered  for  u  and  a  is  sug¬ 
gested  in  Raiffa  and  Schlaifer  (1968)  and  is  called  the  normal-gamma 
density.  This  joint  density  has  the  following  form: 

,  i(v’-l)  -i  0[n'(M-m’)2  +  v'u'] 

h23(u,a^)  0  e  ^  ,  (3.6) 

where  0  =  -»  <  m'  <  u',n',v’  >  0. 

The  joint  prior  density  is  the  product  of  the  densities: 


-  -I-  Bn’  (u-m’  )^ 


h,(u|0)  -  e  ^  —00  <  M  <  - 


and 


-|6v’u’ 

h3(0)  «  e  e "  ,  0  = 


2  ’ 


(3.7) 


(3.7.1) 


where  h2(ljl0)  is  a  normal  probability  density  and  h^(9)  ia  a  gamma 
density.  If  n'  =  0  in  the  joint  density  in  (3.6),  then  p  and  0  have 
independent  priors.  In  this  case,  p  will  have  a  diffuse  prior  and  0 
has  a  gamma'' pr ior .  The  marginal  prior  density  for  p  is  a  studcnt-t  of 
the  form: 


2  -4(v’+l) 

h2(u)  ®  [V  +  (p-m')  n’/u’l 


(3.7.2) 
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Figure  11,  Marginal  posterior  density  for  threshold 
parameter  assuming  beta  prior  density 
for  (J)  over  (50,750)  with  g,  =  2,  g„  =  7. 
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The  advantage  of  using  joint  prior  density  (3,6)  is  that  it  is 
mathematically  tractable.  It  enables  us  to  solve  one  of  the  integrals 
analytically  In  the  evaluation  of  the  posterior  distribution  of  <t>  • 

Let  us  assume  a  prior  density  for  which  is  independent  of  the 
joint  prior  (3.6)  and  denote  it  by  hj^($).  By  using  Bayes'  theorem,  the 
joint  posterior  distribution  of  <)),  U.  and  0  is  obtained  as: 


g(({),U,9  |x) 


*  hj^((i)) 


i=l 


4(n+v’-l) 


-y  fn' (p-m')^+v'u'+  J  [£n(x  -ij!)-ii] 
i=l 


(3.8) 
T 


To  obtain  the  posterior  distribution  of  <J),  (3.8)  is  integrated 
with  respect  to  9  and  p.  The  integral  with  respect  to  9  can  be  evalu¬ 
ated  analytically  and  the  joint  posterior  density  of  ij)  and  p  is  obtained  as 


h^(({)) 


i=l 


^  A  n  , 

y  n'(p-m')  +v'u'+  \  [£n(x.-(}))-p]' 
i=l  ^ 


(3.9) 

y(n'+v'+l) 


The  posterior  distribution  of  ({)  is  obtained  by  integrating  (3.9) 
numerically  with  respect  to  p. 

This  approach  is  applied  to  the  engine  disk  data  by  assuming  a 
beta  prior  density  for  (J). 

A  joint  prior  of  the  form  (3.6)  with  parameters  v'  =  9,  n'  =  10, 

V 

2 

m'  =  4,  u'  =  5  is  assumed  for  o  and  p.  A  beta  prior  density  over 
(50,750)  with  parameters  ~  2  and  =  7  is  assumed  for  <p.  The  pos¬ 
terior  distribution  of  (|i  is  shown  in  Figure  12.  The  mode  of  the 


posterior  distribution  is  at  475  EFH.  Note  that  for  the  independent 

priors  in  Section  3.1,  the  node  was  550  KFH. 

Another  beta  prior  over  (50,750)  with  parameters  =  2  and 

g2  =  5  is  assumed  for  <{).  The  prior  density  is  shown  in  Figure  13. 

2 

The  same  joint  prior  is  assumed  for  y  and  a~ .  By  using  numerical 
integration,  the  posterior  distribution  of  <p  is  obtained  and  is  pre¬ 
sented  in  Figure  14.  The  mode  of  the  posterior  distribution  is  at 

600  EFH.  Note  that  by  assuming  independent  priors  in  Section  3.1  and 

2 

assuming  a  joint  prior  for  y  and  o“  in  this  section,  the  Bayes  estimates 
obtained  for  v>  ^te  similar. 

4.  CONCLL'DI.NC  REMARKS 

In  this  report  we  have  attempted  to  present  a  Bayesian  approach 
for  the  estimation  of  the  threshold  parameter  in  fatigue  crack  models. 
The  nontrivial  forms  of  the  density  functions  have  required  use  of  nu¬ 
merical  techniques  to  evaluate  the  posterior  distributions.  These  nu¬ 
merical  techniques  appear  to  be  efficient  in  evaluating  the  integrals. 

We  have  attempted  to  consider  different  forms  of  priors  in  our  analysis. 
Sometimes  computational  difficulties  have  prevented  us  from  considering 
other  forms  of  prior  densities. 

On  the  basis  of  the  analysis  made,  the  range  (400,550)  EFH  ap¬ 
pears  to  be^a  region  of  high  probability  for  the  threshold  parameter  in 
engine  disks.  Availability  of  more  information  and  selection  of  the 
prior  parameters  appropriately  may  Improve  these  results. 


2 


Figure  14.  Marginal  posterior  density  for  threshold 
parameter  (p  assuming  beta  prior  for  (p 
over  (50,750)  with  gj  =  2,  =  5,  and 

assuming  joint  prior  for  y  and  . 


APPENDIX 


Table  A.l 

Crack.  Initiation  in  10th  Compressor  Disk 


Disk  Number 

Equivalent  Flight  Hours 

1 

1050 

2 

950 

3 

820 

4 

1050 

5 

780 

6 

910 

28 


29 
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